model {
	
	for (l in 1:L){										
		for (j in 1:N) {
			mrgbar[j,l] <- lambda[l] * ( phi[mrg.it[j]]  + theta[mrg.ce[j]] + rho[mrg.c[j]]) 
			mrg.y[j,l] ~ dnorm(mrgbar[j,l], sigma[l])
			}
		}
	

	# Phi 
	for(j in 1:phi.IT){
			phibar[j] <-  phi[phi.z[j]] + 
							  b1[phi.i[j]] * ( -1 + (phi.t[j] - 1) * (2/(phi.T[j]-1)) ) + 
						 	  b2[phi.i[j]] * ( (3*( -1 + (phi.t[j] - 1) * (2/(phi.T[j]-1)) )^2 - 1)/2  ) + 
							  b3[phi.i[j]] * ( (5*( -1 + (phi.t[j] - 1) * (2/(phi.T[j]-1)) )^3 - 3*( -1 + (phi.t[j] - 1) * (2/(phi.T[j]-1)) ))/2  )  
			phi[j] ~ dnorm(phibar[j], sigmainv0[phi.w[j]])
		}

	
	# Lambda / Sigma 
	for(l in 1:L) {
		lambda[l] ~ dnorm(l0, L0)
		sigma[l] ~ dgamma(e0, f0)
		}	
	
	# Theta 
	for(j in 1:mrg.CE){
		theta[j] ~ dnorm(t0, T0)
		}
	
	# Rho
	for(j in 1:mrg.C){
		rho[j] ~ dnorm(r0, R0) 
		}
	
	# Coefficients
	for(j in 1:mrg.I){
			b1[j] ~ dnorm(b0,B0)
			b2[j] ~ dnorm(b0,B0)
			b3[j] ~ dnorm(b0,B0)
		}

}